# clear working environment
rm(list = ls())

#### load the required packages ####
library(tidyverse)
library(PanelMatch)
library(countrycode)

#### load the clean data ####
df <- read_rds("data/dem-transitions-replication-data.rds") %>%
  glimpse()

#### make df into a properly formatted dataframe for PanelMatch ####

# begin coercing df into m
# df needs to be a dataframe class
m <- as.data.frame(df)
# check to see if m is a dataframe
class(m)

# country needs to be an integer to use PanelMatch
m$country <- as.integer(m$country)

# year needs to be an integer to use PanelMatch
m$year <- as.integer(m$year)

# create treatment variable
m$tr <- as.numeric(ifelse(m$v2x_regime >= 2, 1, 0))

# all other variables, except categorical variables, need to be numeric
m$mid <- as.numeric(m$mid)
m$fatal <- as.numeric(m$fatal)
m$ln_cinc <- as.numeric(m$ln_cinc)
m$ln_pecpp <- as.numeric(m$ln_pecpp)
m$terr_dispute <- as.numeric(m$terr_dispute)

# make sure all of the rows are unique
m <- distinct(m, country, year, .keep_all = TRUE)

# create regional democracy indicators
m$un_region <- countrycode(m$country, origin = "cown", destination = "un.region.name")
# match ambiguous historical cases manually
m$un_region <- ifelse(m$country == 260 | m$country == 265 | m$country == 300 | m$country == 315 | m$country == 345 | m$country == 347, "Europe", m$un_region)
m$un_region <- ifelse(m$country == 511, "Africa", m$un_region)
m$un_region <- ifelse(m$country == 678 | m$country == 680 | m$country == 713 | m$country == 730 | m$country == 817, "Asia", m$un_region)

#### effect of democratic transitions on MID participation ####
# create panel data
panel <- PanelData(panel.data = m,
                   unit.id = "country",
                   time.id = "year",
                   treatment = "tr",
                   outcome = "mid")

# create matched sets
PMresults.CBPSweight.5.5.mid.att <- PanelMatch(panel.data = panel, lag = 5, refinement.method = "CBPS.weight",
                                               match.missing = TRUE, exact.match.variables = c("un_region"),
                                               covs.formula = ~ I(lag(mid, 1:5)) + I(lag(ln_cinc, 1:5)) + I(lag(ln_pecpp, 1:5)) + I(lag(terr_dispute, 1:5)),
                                               qoi = "att", forbid.treatment.reversal = FALSE,
                                               use.diagonal.variance.matrix = TRUE)

# look at the pretreatment balance
pre1 <- get_covariate_balance(panel.data = panel,
                              PMresults.CBPSweight.5.5.mid.att,
                              covariates = c("mid"))

# format for plotting
pre1 <- as.data.frame(pre1) %>%
  rownames_to_column(var = "year") %>%
  mutate(year = str_replace(year, "_", "-")) %>% # this replaces an underscore with a minus sign
  mutate(model = "All MIDs") %>%
  filter(year != ("t-0")) %>% # excludes the year the transition occurs because of timing within year uncertainty
  glimpse()

#### effect of democratic transitions on fatal MID participation ####
# create panel data
panel <- PanelData(panel.data = m,
                   unit.id = "country",
                   time.id = "year",
                   treatment = "tr",
                   outcome = "fatal")

# create matched sets
PMresults.CBPSweight.5.5.fatal.att <- PanelMatch(panel.data = panel, lag = 5, refinement.method = "CBPS.weight",
                                                match.missing = TRUE, exact.match.variables = c("un_region"),
                                                covs.formula = ~ I(lag(fatal, 1:5)) + I(lag(ln_cinc, 1:5)) + I(lag(ln_pecpp, 1:5)) + I(lag(terr_dispute, 1:5)),
                                                qoi = "att", forbid.treatment.reversal = FALSE,
                                                use.diagonal.variance.matrix = TRUE)

# look at the pretreatment balance
pre2 <- get_covariate_balance(panel.data = panel,
                              PMresults.CBPSweight.5.5.fatal.att,
                              covariates = c("fatal"))

# format for plotting
pre2 <- as.data.frame(pre2) %>%
  rownames_to_column(var = "year") %>%
  mutate(year = str_replace(year, "_", "-")) %>% # this replaces an underscore with a minus sign
  mutate(model = "Fatal MIDs") %>%
  rename(mid = fatal) %>%
  filter(year != ("t-0")) %>% # excludes the year the transition occurs because of timing within year uncertainty
  glimpse()

pre <- rbind(pre1, pre2)

#### create parallel trends plot for Figure F.1 Panel A ####
id_panel_a <- ggplot(data = pre, aes(x = reorder(year, desc(year)), y = mid)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point() +
  geom_line(aes(group = 1)) +
  facet_wrap(~model) +
  labs(y = "Standardized Difference in Pr(Conflict)", 
       x = "Years before Democratic Transition",
       title = "Panel A: Descriptive Differences in Pr(Conflict) before Transition") +
  ylim(-2, 2) +
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(color = "black", size = 10),
        axis.title.y = element_text(size = 10),
        legend.title = element_text(size = 10))


#### save ggplot object to combine with placebo test later ####
write_rds(id_panel_a, "data/identification-assumptions-panel-a.rds")
